function result = miecoated_ab3(m1, m2, x, y)

  % Computation of Mie Coefficients, a_n, b_n,
  % of orders n=1 to nmax, complex refractive index m=m'+im",
  % and size parameters x=k0*a, y=k0*b where k0= wave number
  % in the ambient medium for coated spheres, a=inner radius,
  % b=outer radius m1,m2=inner,outer refractive index;
  % p. 183 in Bohren and Huffman (1983) BEWI:TDD122
  % but using the bottom equation on p. 483 for chi_prime
  % C. Mätzler, August 2002

  m = m2 / m1
  nmax = round(2 + y + 4 * y.^(1/3))
  n = (1:nmax)
  nu = (n + 0.5)
  u = m1 * x
  v = m2 * x
  w = m2 * y
  su = sqrt(0.5 * pi * u)
  sv = sqrt(0.5 * pi * v)
  sw = sqrt(0.5 * pi * w)
  sy = sqrt(0.5 * pi * y)
  pu = su .* besselj(nu, u)
  py = sy .* besselj(nu, y)
  pv = sv .* besselj(nu, v)
  pw = sw .* besselj(nu, w)
  p1u = [sin(u), pu(1:nmax - 1)]
  p1y = [sin(y), py(1:nmax - 1)]
  p1v = [sin(v), pv(1:nmax - 1)]
  p1w = [sin(w), pw(1:nmax - 1)]
  ppv = p1v - n .* pv ./ v
  ppw = p1w - n .* pw ./ w
  ppy = p1y - n .* py ./ y
  chv = -sv .* bessely(nu, v)
  chw = -sw .* bessely(nu, w)
  chy = -sy .* bessely(nu, y)
  ch1y = [cos(y), chy(1:nmax - 1)]
  gsy = py - i * chy
  gs1y = p1y - i * ch1y
  gspy = gs1y - n .* gsy ./ y
  du = p1u ./ pu - n ./ u
  dv = p1v ./ pv - n ./ v
  dw = p1w ./ pw - n ./ w
  chpw = chw .* dw - 1 ./ pw
  uu = m .* du - dv
  vv = du ./ m - dv
  pvi = 1 ./ pv
  aaa = pv .* uu ./ (chv .* uu + pvi)
  bbb = pv .* vv ./ (chv .* vv + pvi)
  aa1 = ppw - aaa .* chpw
  aa2 = pw - aaa .* chw
  bb1 = ppw - bbb .* chpw
  bb2 = pw - bbb .* chw
  aa = (py .* aa1 - m2 .* ppy .* aa2) ./ (gsy .* aa1 - m2 .* gspy .* aa2)
  bb = (m2 .* py .* bb1 - ppy .* bb2) ./ (m2 .* gsy .* bb1 - gspy .* bb2)

  result=[aa; bb]
